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ABSTRACT 

Atmospheric refraction affects to various degrees exoplanet transit, lunar eclipse, 
as well as stellar occultation observations. Exoplanet retrieval algorithms often use 
analytical expressions for the column abundance along a ray traversing the atmosphere 
as well as for the deflection of that ray, which are first order approximations valid for 
low densities in a spherically symmetric homogeneous isothermal atmosphere. We 
derive new analytical formulae for both of these quantities, which are valid for higher 
densities, and use them to refine and validate a new ray tracing algorithm which can be 
used for arbitrary atmospheric temperature-pressure profiles. We illustrate with simple 
isothermal atmospheric profiles the consequences of our model for different planets: 
temperate Earth-like and Jovian-like planets, as well as HD189733b, and GJ1214b. 
We find that, for both hot exoplanets, our treatment of refraction does not make 
much of a difference to pressures as high as 10 atmosphere, but that it is important 
to consider the variation of gravity with altitude for GJ1214b. However, we find that 
the temperate atmospheres have an apparent scale height significantly smaller than 
their actual density scale height at densities larger than 1 amagat, thus increasing 
the difficulty of detecting spectral features originating in these regions. These denser 
atmospheric regions form a refractive boundary layer where column abundances and 
ray deflection increases dramatically with decreasing impact parameter. This refractive 
boundary layer mimics a surface, and none of the techniques mentioned above can 
probe atmospheric regions denser than about 4 amagat on these temperate planets. 

Key words: atmospheric effects - methods: analytical methods: numerical planets 
and satellites: atmospheres - radiative transfer. 


1 INTRODUCTION 

Most analytical expressions of an exoplanetary primary 
transit lightcurve are derived from analytical expressions of 
the optical depth along an unrefracted ray, i.e. that travels 
in a straight light, through a homogeneous isothermal at¬ 
mospheres (see e. g. Fortney 2005; Lecavelier des Etangs et 
al. 2008; Benneke & Seager 2012; Howe & Burrows 2012; 
de Wit & Seager 2013; Griffith 2014). However, several pa¬ 
pers (Hui & Seager 2002; Sidis & Sari 2010; Garcia Munoz 
et al. 2012; Betremieux & Kaltenegger 2013, 2014; Misra, 
Meadows & Crisp 2014; Misra & Meadows 2014) have em¬ 
phasised the importance of the refractive bending of light 
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by exoplanetary atmospheres on these lightcurves. Refrac¬ 
tion is also important in various observational geometries in 
our Solar System. It influences the perceived brightness of a 
star during a stellar occultation by a planetary atmosphere 
(see the review by Smith & Hunten 1990) which can then be 
used to infer the planet’s atmospheric composition. It also 
determines the brightness of moons eclipsed by their parent 
planet. Indeed, as a moon moves deeper in the penumbra, 
solar radiation must traverse deeper atmospheric regions of 
the eclipsing planet in order to be bent sufficiently to reach 
the moon. Several Lunar eclipse observations were analysed 
to derive the transmission of Earth’s atmosphere (Palle et 
al. 2009; Vidal-Madjar et al. 2010; Garcia Munoz et al. 2012, 
Arnold et al. 2014) using this principle. 

Given the importance of atmospheric refraction in these 
various observational geometries, it is worthwhile to have an- 
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alytical formulae which include refraction, both to compute 
the ray’s deflection, as well as the integrated number density 
along the deflected ray, or its column abundance. Combined 
with an atmosphere-averaged extinction cross section, one 
can compute the optical depth along the ray, and the cor¬ 
responding atmospheric transmission with Beer’s law. How¬ 
ever, computing the column abundance and the deflection 
requires solving complicated integrals for which analytical 
expressions exist only under simplifying assumptions. In sec¬ 
tion [2] we derive analytical solutions of these integrals by 
expanding the integrand in a Taylor-series up to the sec¬ 
ond order with respect to refractivity, and then evaluate the 
integral for each of these terms. This is done for a homoge¬ 
neous isothermal atmosphere, since this is the level of com¬ 
plexity considered by some exoplanet atmosphere retrieval 
algorithm (see e.g. Charbonneau et al. 2009; Howe & Bur¬ 
rows 2012; Anglada-Escude et al. 2013; Ehrenreich et al. 
2014; Waldmann et al. 2014). We also derive an order-of- 
magnitude estimate of the sum of the remaining uncom¬ 
puted higher-order terms in the Taylor series, in order to 
determine the density region over which our analytical ex¬ 
pressions are useful. 

Since planetary atmospheres are usually not isothermal, 
we have to rely on numerical methods to deal with these 
general cases for which analytical solutions do not exist. In 
section [3j we describe MAKEXOSHELL, a new numerical 
ray tracing algorithm that computes the column abundance 
and the deflection for rays traversing a spherically symmet¬ 
ric atmosphere given an arbitrary one-dimensional (1-D) 
temperature-pressure profile as input. We also describe how 
we compute the model atmosphere, as well as its refractiv¬ 
ity, which are inputs to MAKEXOSHELL. In section |4l we 
compare the output of the numerical algorithm with our an¬ 
alytical expressions for an Earth-like planet with a homoge¬ 
neous isothermal atmosphere, in order to mutually validate 
both methods. We also discuss the impact of our new model 
on temperate Earth-like and Jovian-like exoplanets, as well 
as for two already well-observed exoplanets: HD189733b and 
GJ1214b. We summarise our results in section 0 


2 ANALYTICAL EXPRESSIONS 
2.1 Basic equations 

As light rays traverse a planetary atmosphere, they are bent 
by refraction toward the surface of the planet (see Fig. |T]) 
due to the exponential increase of atmospheric refractivity 
(the refraction index minus one), with decreasing altitude. 
The trajectory of a ray is described by an invariant, b, which 
is equal to its impact parameter (Phiimey & Anderson 1968; 
Betremieux & Kaltenegger 2013, 2014) 

b = r(l + v) sin 6 = ro(l + vq) = Rtop sin#o. (1) 

Here, r is the radial position with respect to the centre of the 
planet (i.e. radius + altitude), which for the sake of brevity 
we will refer to as height throughout the paper, while v and 9 
are the refractivity of the atmosphere and the zenith angle 
of the ray at that height, respectively, ro is the minimum 
height reached by a ray, or its grazing height, where the 
atmospheric refractivity is vo. Rt. op is the radial position of 
the top of the atmosphere and 9o is the zenith angle of the 
ray at the top of the atmosphere, or simply its incidence 


TO OBSERVER A(fl TO STAR 



Figure 1. Trajectory through a planetary atmosphere of an ob¬ 
served light ray during an exoplanetary primary transit. The solid 
body, or opaque region, of the planet (dark grey) and its atmo¬ 
sphere (light grey) are shown, along with the axis running from 
the observer to the centre of the planet (dashed line) with the 
observer to the left and the star to the right. The radius of the 
planetary surface, R p , and of the top of the atmosphere, Rtop, as 
well as the zenith angle, 9, of the ray for a given height, r, are 
also indicated. A light ray observed with an impact parameter b, 
entered the atmosphere with an incidence angle 9q a projected 
distance b' to the centre of the planet with respect to the ob¬ 
server, reached a grazing height ro, and is deflected by uj by the 
atmosphere. 


angle. The atmosphere deflected the ray by oo, and the ray 
entered the atmosphere a projected distance b' to the centre 
of the planet, with respect to the observer, given by 


b' = ( b . 
\ sin 9 o ) 


( 2 ) 


The change in zenith angle of the ray with height is 
given by, 


d9 tan 9 f dv\ tan 9 
dr (1 + u) \ dr ) r 


(3) 


The first term is the refractive term and comes from Snell’s 
law. The second term is the geometric term, and comes from 
the equation of an unrefracted ray, b = rsinfl, which pro¬ 
ceeds in a straight line. Note that all of the equations involv¬ 
ing angles are valid when these are expressed in radians, even 
though in section [4] we quote results in degrees. 

The atmospheric refractivity is given by, 

V = ( H y: fjVSTPj = ( n ^ vs TP, (4) 

V nSTP ) V n STP J 


where n is the number density, tlstp is the number density 
at standard temperature and pressure (STP) also known as 
Loschmidt’s number, fj is the mole fraction of the j th chem¬ 
ical species, vsrpj is the STP refractivity of the j th species, 
while vstp is the STP refractivity of the atmosphere. The 
ratio inside the parenthesis in equation[I]is the number den¬ 
sity expressed in units of amagat. Note that the STP refrac¬ 
tivity is simply the refractivity for a density of 1 amagat, 
which can be achieved for temperatures other than 273.15 K 
with the appropriate pressure. 

The number density, n, in an atmospheric region varies 
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with height, and this dependence can be expressed by 


nu =n L e- (ru ~ rL)/H , 


(5) 


so that 


dn n 
dr = ~H' 


( 6 ) 


where H is the density scale height in that atmospheric re¬ 
gion, and the U and L subscripts denote quantities for the 
upper and lower boundaries of the region, respectively. Thus, 


dv 

vstp , 

7dn\ 

f n. \ 

i vstp 

V 

dr 

nsTP ' 

\dr) 

V nsTP) 

1 H 

H 


(7) 


Note that this result holds only for a homogeneous atmo¬ 
spheric region where the composition of the major con¬ 
stituents that contribute significantly to the refraction in¬ 
dex does not change with height, i.e. that vstp is constant 
with height. This assumption will be carried throughout all 
derivations. 


2.2 Lower atmospheric boundary 


In order for a ray to traverse the atmosphere, 9 must increase 
as r decreases (or dd/dr < 0) so that a downward-going ray 
reaches a minimum height, ro, where it grazes the atmo¬ 
sphere (9 = 7t/ 2). If dd/dr > 0, the ray will spiral deeper 
into the atmosphere until it is absorbed. Using equations [3j 
we can derive that 


r (dv\ 

(1 + v) V*v < 


( 8 ) 


is the necessary condition for a ray to traverse the atmo¬ 
sphere. Using equation [7] this can be rewritten as, 


1 + v 


r 

H <1 - 


(9) 


In the special case when this expression is equal to one, a ray 
which is tangential to the atmosphere will circle the planet 
as its radius of curvature matches its height. Knowing the 
height-dependence of the density, and thus refractivity, this 
equation can be solved recursively for this special height. As 
a ray approaches this height, both the column abundance 
and the deflection tend toward infinity. Atmospheric regions 
below this height cannot be probed by transmitted stellar 
radiation under any circumstance during a stellar occulta- 
tion or an exoplanet primary transit. This is true only for 
a perfectly spherical planet free of horizontal density gradi¬ 
ents (i.e. no weather or turbulence). Departure from these 
assumptions creates a diffuse boundary about which a graz¬ 
ing ray can or cannot escape the planet depending on local 
conditions. Furthermore, since refractivities vary with wave¬ 
length, so does the height of this lower boundary. 

This lower boundar y is always located b elow th e crit i- 
cal altitude, described in lBetremieux fe Kaltenegged (120141 s ) . 
below which atmospheric regions cannot be probed by trans¬ 
mission spectroscopy when an exoplanet occults the central 
area of its host star. This critical altitude occurs where a ray, 
that comes from the opposite limb of the star, reaches an 
atmospheric critical density which causes it to be deflected 
by a critical deflection, u> c , toward the observer. During the 
course of the transit, the symmetry is broken and part of 
the planetary limb can be probed to lower altitudes while 
the opposite limb’s observable region is restricted to higher 


altitudes. However, the lowest altitude that one can probe 
only asymptotically approaches the lower boundary deter¬ 
mined by equation [51 so that this boundary is the deepest 
atmospheric region that can ever be probed by exoplanet 
transmission spectroscopy or by stellar occultation. 


2.3 Column abundance and ray deflection 
integrals 

The column abundance, N, along a refracted ray, as well as 
its deflection, u>, can be computed from 

POO POO 

N = 2 dN = 2 F\dr (10) 

•>r o J r 0 

and 

POO POO 

oj = 2 dcu = 2 F2dr ( 11 ) 

Jro Jr 0 

where the factor of 2 comes from a ray traversing an altitude 
region twice. Note that we have labelled the integrands in 
terms of the integration parameter r in both equations for 
future use (see section E3J. 

The incremental column abundance, dN, of a ray is 
given by, 

T) 

dN = - -dr. (12) 

cos 9 

The incremental deflection of a ray (expressed in radians), 
dw, is simply the first term in equation[3](see also Goldsmith 
1963; Auer & Standish 2000), which we rewrite as 


did = — 


tan 9 

(T+^y 



u 

1 + v 


tan 9 
H 


dr 


(13) 


for a homogeneous planetary atmosphere. 

The trigonometric functions can be rewritten entirely 
in terms of r and v. From equation [Tj we have 


. b r 0 (l + v 0 ) 

sm 9 = —- = —)- 

r(l + v) r( 1 + v) 


(14) 


1 = 1 = (\_( r o( 1 + i/ o) ^ 2 \ 1 

cos 9 y'l _ sin 2 0 ^ \r{l + v) ) ) 

and 


tan 9 


sin 9 
cos 9 


( f r(l + v) y 
\\ro(l + vo)J 



(15) 


(16) 


Equations m and ED are then given, for a homogeneous 
atmosphere, by 



where most of the height-dependence is buried inside the 
refractivity. 
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2.4 Taylor-series expansion of integrands 

No analytical solutions exist for the complicated integrals 
in equations El and El even in the case of a homogeneous 
isothermal atmosphere. However, refractivities are typically 
small, on the order of ICC 4 at STP, such that one can expand 
the integrand in both equations in Taylor series in uo, prior 
to carrying-out the integration. We carry out this Taylor 
expansion for both integrands F\ and F 2 to the second order 
such that 


dN = dN 0 + dm + dN 2 + ... 

dFY 


= F 1 U =0+V0 [ — 


"0=0 




and 

du 


■ T dcu 1 T duj2 


d?F\ 

dvl 


dr F 2 
dv 0 


i / 0 =0 


(19) 

+ 

( 20 ) 
+ .... 


i'0=0 


For a homogeneous isothermal atmosphere, we have the 
additional relations, 


— (r — rn) / H 

n = noe v U7/ 


( 21 ) 


and 

—(r— rn)/H 

v — uoe v U7/ 


( 22 ) 


where no and uq are the number density and the refractivity 
at a ray’s grazing height, respectively. The derivation of the 
various Taylor expansion terms is rather tedious and lengthy, 
so we will here only quote the results: 


Fi|i/ 0 =o= 


n _ n(r/ro) 
\/l- (r 0 /r) 2 V( r /r o) 2 - 1 


(23) 


dF\\ 


dv °J „ 0 =0 


n(r/ro)(l-e- (r - ro)/H ) 
{{r/r o) 2 - I ] 3 / 2 


(24) 


3n(r/r 0 ) 3 (l - e 


3(1 _ „-(>— r 0 )/H\2 


du n / „ 

U / 1/0=0 


and, 

F2\v o =0= 0 


[(r/r 0 ) 2 — I ] 5 / 2 

2n(r/r 0 )(l-e- {r - ro)/H ) 
[(r/r 0 ) 2 - I ] 3 / 2 


(25) 


(26) 


-(’— r o)/ H 


dF 2 \ _ 1 

du ° ) i/ 0 =o H \Z{r/ro) 2 - 1 


(27) 


d 2 F 2 

di/ 2 


a -(r-r 0 )/H 


H 


^0=0 “ V ( r / r o ) 2 ~ 1 

(r/r 0 ) 2 (l-e- < ~ r - ro)/H ) 


o -(r-r 0 )/H 


(r/r 0 ) 2 - 1 


(28) 


2.5 “Modified Bessel functions” solution 

At first glance, it seems that we have made the problem 
more complicated as we have transformed each of our initial 
integrals into sums of three integrals, each of which also 


seem impossible to solve. However, a suggestion of a solution 
appears when on e cons id ers t he modified Bessel functions of 
the second kind (|Arfkenlll985i 'l. 

*■<*> - (TTW (2) l < 2! » 

of which the first two functions are 

/ OO 

e~ zx (x 2 — l)~ 1 ^ 2 da:, (30) 

and 


/ OO 

e~ zx (x 2 — l) 1 ^ 2 dx, (31) 

which were derived from equation 1291 knowing that 
(-1/2)!= tt 1/2 . 

Consider, as a first step, the zeroth term in the column 
abundance integral, 


No 



dNo = 2 



n(r/r 0 ) 

\](r/ro) 2 - 1 


dr 


no{r/ro)e- (r - ro)/H 
\J{r/r o) 2 - 1 


dr. 


(32) 


With the change of variable, x = r/ro, it becomes 


f°° re ~{r 0 /H)x 

No = 2noroe o/ / _ 1 ) 1/2 dx ( 33 ) 

which can be integrated by parts to yield 


No = [2n 0 roe r ° /Jf (a: 2 - l) 1 / 2 e -( r * , / H )“] / 

/ OO 

e ~ < . r o/ H ') x ^ x 2 _ i) 1 ' /2 da:. (34) 

The values of the first term at 1 and 00 are both zero so 
the first term drops out. Whereas it is straightforward to 
derive for the value of 1, the 00 case requires the use of 
l’Hopital’s rule. The second term can be rewritten in terms 
of the modified Bessel functions of the second kind, so that 
we obtain 

No = 2noroe a/H K 1 (r 0 /H) = 2 n 0 r 0 Kt ( r 0 /H) (35) 

where we adopt the following short-hand notation, 

K;(y) = e y K t (y) (36) 


for the remainder of this paper. This expression for the col¬ 
umn abundance has been previously derived in terms of 
these modified Bessel functions (see e.g. Griffith 2014). How¬ 
ever, this is only the leading term in our Taylor series, and 
corresponds to the non-refractive case. 

Each of the terms in the column abundance and de¬ 
flection expansion can also be rewritten in terms of these 
modified Bessel functions using the same change of variable, 
and then integrating by parts a number of times. At each 
integration by parts, we get a term that drops out when eval¬ 
uated at 1 and 00 and an integral term that can either be 
expressed as a modified Bessel function or must be further 
integrated by parts. Again, the derivation is straightforward 
but rather lengthy so we only here quote the results: 


m = 2n 0 r 0 


m 



(37) 
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—“(t) * T T 


-T«(f) + 4*(2.) + *(g) 


* ( 3 r 0 


2 r 0 


and 
cjo = 0 


(38) 

(39) 



Note that the leading term of the Taylor-series expansion for 
deflection, shown in equation [31 is 0 as expected since the 
leading term is the non-refractive case. Without refraction, 
light is not bent by the atmosphere. 


2.6 “General power-series” solution 

Expressions in po wer-se ries exists for these modified Bessel 
functions ( Arfkcnl ll985i h and are given by: 


Kt{y) 



(4t 2 - 1) 
8 y 


(At 2 - l)(4t 2 

128t/ 2 



In our case, y is always a multiple of ro/H, which is al¬ 
ways large. Thus, the second, third and higher-order terms 
in equation El decrease by successive powers of multiples of 
ro/H, and are significantly smaller than the leading term in 
this series. 

Collecting terms, equations 1191 and 1201 can be rewritten 
more generally as 


N = 

and 

oj = 


2-irro 

H 


(■ n 0 H) 


OO 

i + a + E^tf)' 


3-1 


2nro 


H 


-vo 


i + n> + £>(tt)' 


3 =1 


(43) 


(44) 


where Cj and Dj are coefficients for the column abundance 
and deflection, respectively. Even though we have only ex¬ 
plicitly collected terms up to j = 2 for the column abun¬ 
dance, and j = 1 for the deflection, we can determine by 
inspection that higher-order terms follow this trend. These 
coefficients are themselves power-series expansion in term of 
H/ro , but each of them quickly converges since this ratio is 
always small. 

Both series depends on 2-Kro / H, the ratio of gas col¬ 
umn abundance along a grazing ray to that along a ver¬ 
tical ray (given by noH) for a curved exponential atmo¬ 
sphere. This ratio, which we will call the slant factor, has 
been described before (see e.g. Fortney 2005). The leading 
non-zero term in both series matches previous expressions 
of the column abundance without refraction (see references 


listed in section [l]), as well as first-order expressions for the 
ray deflection (Goldsmith 1963), respectively. Higher-order 
terms in the series are surprisingly dependent on powers of 
( roVo/H ), and not simply v o- Whereas uo is usually very 
small, it increases exponentially with depths so that com¬ 
bined with ro/H , which is large, ( roVo/H ) can become of 
order unity. Hence, higher-order terms in the series are far 
from negligible when a ray travels deep enough in a plane¬ 
tary atmosphere. So over what range of densities can we use 
our analytical expressions? 

We have only computed the coefficients up to the sec¬ 
ond order in refractivity as the length of the calculations, 
and correspondingly the odd of making a mathematical mis¬ 
take, increases steeply with order. However, we can make 
an order-of-magnitude estimate of the sum of the terms for 
which we have not explicitly computed the Cj and Dj coef¬ 
ficients, for both the column abundance, N res t, and the de¬ 
flection, LJrest- This is done by assuming that the unknown 
coefficients are all identical to the highest-order computed 
coefficient. Under this crude assumption, we can write 

f^)' (*) 

1=0 

and 

/27rro frov 0 \ 2 ^(r 0 vo\ l , . 

UJrest - y —^~V 0 D 1 [-fi-J {-ff~) ( 46 ) 

1=0 

where the leading term for C 2 and D\ are 0.27 and 0.41, 
respectively. The infinite sum, which appears in both equa¬ 
tions EH and sa is a geometric series and converges to 


E 


/ r- 0 zx 0 \ 

1 C r 0 v 0 \ 

\ H ) 

~ V H J 


(47) 


but only when 


rpvp 

H 


1. 


(48) 


It is interesting to note that the condition for the conver¬ 
gence of this series is very close to the condition for a ray 
to traverse the atmosphere, shown in equation E] Indeed, 
they only differ by a factor (1 + vq) when evaluated at ro. 
Since (1 + vo) > 1, equation[9]is automatically satisfied when 
equation EH1 is satisfied. The fact that they are not exactly 
the same is probably due to our crude method of estimating 
the remaining terms. 


3 NUMERICAL MODEL 
3.1 Ray-tracing algorithm 

We have made substantial improvements to an algorithm 
that, given a 1-D model atmospheric temperature-pressure 
altitude profile, traces rays through that atmosphere and 
computes the number density column abundance along the 
rays as well as their deflection by the atmosphere. This algo¬ 
rithm, which has been previously described by Betremieux 
& Kaltenegger (2013, 2014) and references therein, suffered 
from numerical instabilities in the computed ray deflection, 
which started at densities of about 5 x 10 -5 amagat, and be¬ 
came worse with lower densities, such that the algorithm was 
useless for densities below 3 x 10 -6 amagat. These numeri¬ 
cal instabilities were not only due to improper treatment of 
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integer to double precision conversion of variables, but also 
intrinsically due to the method. Indeed, we were computing 
the ray deflection by first computing the azimuthal travel of 
the ray through the atmosphere, A (j> (see Fig. [TJ, and then 
using 

w = A0 + 200 - 7T, (49) 

thus subtracting two large numbers from each other to get a 
much smaller quantity. We have modified the algorithm and 
reduced these numerical instabilities considerably. In this 
section, we will describe the new algorithm which we name 
MAKEXOSHELL. 

The inputs to MAKEXOSHELL include a 1-D model 
atmosphere’s temperature and pressure, or T and P, as a 
function of altitude, z. Also required are the radius of the 
planet, Rp, the highest and lowest altitudes of the compu¬ 
tational region, the number of rays and layers sampling this 
altitude region, and the STP refractivity of the atmosphere. 
Given temperatures and pressures from the model atmo¬ 
sphere, typically in increments of 1 km altitude, MAKEX¬ 
OSHELL uses the ideal gas law to determine the correspond¬ 
ing densities. This serves as a basis for the rest of the compu¬ 
tation done on a much finer altitude grid, or computational 
grid. 

MAKEXOSHELL assumes that the atmosphere within 
a model atmosphere layer follows an exponential behaviour 
with a constant density scale height as displayed in equa¬ 
tion [5] Once the scale height within that layer is computed, 
using the densities at its upper and lower boundaries, densi¬ 
ties are interpolated on the computational grid within that 
layer. Thus, MAKEXOSHELL builds-up a density profile 
on the computational grid from the model atmosphere den¬ 
sities. This is then combined with the input vstp, as per 
equation [4] to compute the altitude-dependent refractivity 
of the atmosphere in the middle of each computational layer. 

Rays are spread uniformly in altitude, in terms of their 
grazing height, across the specified computational region. At 
each of these grazing heights, the density, no, and refractiv¬ 
ity, vo are determined from which the impact parameter of 
each ray can be computed according to equation[l] With the 
impact parameter of a ray, b, as well as the density and the 
refractivity of the atmosphere on the computational grid, it 
is possible to compute the column abundance of gas inter¬ 
cepted by a ray, as well as the ray’s deflection through the 
atmosphere. 

The computational grid is made up of layers that are 
thin enough that their densities, and hence their refractivi- 
ties, are assumed constant across each layer. Thus, density 
and refractivity changes occur at their boundaries. In such 
a case, rays follow a straight line within each computational 
layer, as shown in Fig. [2] We can see from this simple ge¬ 
ometry that 

d(f> = 9 l — 9jj, (50) 

and 

L = ru sin d(j> = ds sin 9l, (51) 

where dcj> and ds are the azimuthal travel of the ray and the 
path length of the ray within the computational layer, re¬ 
spectively. Subscript U and L denote quantities evaluated at 
the upper and lower boundary of the computational layer, 
respectively. Quantities without a subscript are evaluated 



Figure 2. Geometry of a ray travelling through an atmospheric 
layer with constant refractivity, i.e. in a straight line. 


in the middle of the computational layer, except for the 
impact parameter and the scale height. The former is ray- 
dependent, while the latter is determined by the model at¬ 
mosphere layer to which the computational layer belongs. 
Both ds, given by 

j „ _ r u s in(d0) r v sin(0 L - 0 V ) /crrA 

• /) ft ' 

sin ul sin0L 

derived from equation 15II and d(f> depend on the zenith an¬ 
gles of the upper, 9u, and lower boundaries, 9l, of the com¬ 
putational layer. These are computed with 


sin 0(7 = 


b 

(1 + v)r v ’ 


(53) 


and 


sin0 L = b —, (54) 

(1 + v)r L 

which have been derived from equation [l] 

The column abundance intercepted by a ray in a com¬ 
putational layer is derived simply by 


dN = nds 


(55) 


while its deflection is given by 

duJ= ih (i ) d4> • (56) 

Equation [SS] is derived from equation 1131 recognising that 

71 tal1 ^ , 

d(j> = - dr, (57) 

r 

which is the geometric term in equation[3] MAKEXOSHELL 
integrates ds and d(f> across all computational layers tra¬ 
versed by each rays, in order to get the column abundance 
and the deflection across an upward or downward pass for 
each rays. The total column abundance and deflection ex¬ 
perienced by a ray is twice that quantity since the rays goes 
through these layers twice: once moving downward, and once 
moving upward through the atmosphere. 

Although the computational grid used in Betremieux & 
Kaltenegger (2013, 2014) was composed of 0.1 km-thick lay¬ 
ers, we have found that the precision of the computations 
improves significantly when the thickness of the computa¬ 
tional layers decreases (see section iTTll . and that the desired 
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thickness depends on the atmospheric scale height. Hence, 
we now adjust this parameter as required by the atmosphere 
for which we do the computations. 

The contributions to the column abundance, dN top , and 
ray deflection, dujtop, from atmospheric layers above the top 
of our computational altitude region are computed by as¬ 
suming that these upper atmospheric regions are compressed 
into a shell one scale height thick. At these low densities, the 
ray trajectory is essentially straight, so that the path length 
of the ray, dstop, inside this shell, is given by 

dstop = \J Rtop cos 2 0' Q + 2H top Rtop + H? op -R t0 p cos 0' o (58) 

where the subscript top refers to quantities evaluated at the 
top of our computational region, or top of the atmosphere, 
with a corresponding refractivity, vtop, and zenith angle, 9' 0 
(see also Fig. [T]). Then, the contribution from these upper 
layers to the column abundance is simply given by 

dNtop — nt.opdstop, (59) 

and to the deflection of the ray by 

duhop = ftop ( \ sin d' 0 . (60) 

\Jdtop J 

The latter equation was derived by using the small angle 
approximation in equation M to derive drf>, which is then 
substituted into equation Eg 


3.2 Refractivities 


One of the critical parameters for computing the deflec¬ 
tion of the ray as well as the column abundance along the 
refracted ray is the STP refractivity of the atmosphere, 
which depends on the STP refractivity of its major chemi¬ 
cal species (see equation [4]). In this paper, we use the same 
refracti v ities for_ N 2 , O 2 , C O 2 , and Ar, as in Table 2 from 
iBetremieux fc Kalteneggeil (l2013h . Note that t his t able has 
a small error: the STP refractivity for O 2 from lBates! (Il984l l 
can be computed for wavelengths above 0.546 pm, contrary 
to what this table states. We are also considering two more 
species: H 2 and He. The STP refr activity of Helium i s com - 
puted with the expressio n from iMansfield fe Peckl d 196914 . 
described in i Weben (l2003l l . For molecular hydrogen, we start 
from the formulation by Ford & Brown (1973) of the H 2 
Rayleigh cross section as a function of its rotational quan¬ 
tum number J, from which we compute the cross section, 
an, for an equilibrium H 2 population at 300 K. We then 
determine its STP refractivity by using Rayleigh’s formula, 


vstp = 


nsTP / 3 or 
W 2 V 327T 3 ’ 


(61) 


where w is the wavenumber of the radiation. 

Since the refractivity of molecules hardly changes in the 
near-infrared (NIR) and infrared (IR) part of the spectrum, 
we can pick a single refractivity to represent that entire 
spectral region with minimal errors. In this paper, we use 
the refractivity of various molecules at 1 pm (shown in Ta¬ 
bic [J ). Note that th e deriv ation of the refractivity of H 2 
fron uFord fe Brownel (Il973i 'l includes a very slight tempera¬ 
ture dependence (about a 0.3 per cent change in refractivity 
from 273.15 to 2000 K), but we will ignore it here and use 
the refractivity in Table [I] 


Table 1. STP refractivities of various gases and planetary atmo¬ 
spheres at 1 pm (see sect.ion If 21 for references). 


Gas / Atmosphere 

VSTP 

h 2 

1.37 x icr 4 

He 

3.48 x icr 5 

n 2 

2.95 x icr 4 

O 2 

2.68 x icr 4 

Ar 

2.79 x icr 4 

CO 2 

4.43 x icr 4 

Jupiter 

1.23 x icr 4 

Earth 

2.90 x 10 -4 


We can compute the STP atmosphere refractivity for an 
Earth-like atmosphere, considering the contribution of N 2 , 
O 2 , Ar, a n d CO 2 , using the a tmosphere composition listed in 
iLodders fe Feelev, Jr.l Jl998T ). For a Jovian-like atmosphere, 
we use the mole fraction of he lium (0.1357) measured b y th e 
Galileo probe, also listed in ILodders fe Feglev. Jr.l (11998T) . 
and assume that the remainder of the atmosphere is com¬ 
posed of molecular hydrogen. 


3.3 Atmospheric models 


MAKEXOSHELL requires as input an atmospheric tem¬ 
perature and pressure profile as a function of altitude. In 
this paper, we only consider homogeneous isothermal atmo¬ 
spheres, so that both the temperature, T, and the compo¬ 
sition are constant with altitude, which is currently done 
by some exoplanet atmosphere retrieval algorithm (see e.g. 
Charbonneau et al. 2009; Howe & Burrows 2012; Anglada- 
Escude et al. 2013; Ehrenreich et al. 2014; Waldmann et al. 
2014). 

To build the model atmosphere, we start by specifying 
the mass of the planet, Mp, the surface pressure, P a , the ref¬ 
erence pressure, Pp, at the planetary radius, Rp, the thick¬ 
ness of the atmosphere, AZ a tm, and the isothermal tempera¬ 
ture, T, of the atmosphere. We compute the mean molecular 
mass, m, from our specified bulk atmospheric composition. 
We discretize our atmosphere in small altitude increments, 
Az. We then determine the pressure scale height, H[, at each 
altitude by 


H' = 


kT 

mgt 


(62) 


where the gravitational acceleration, gi, can either vary with 
altitude 


9i = 9 s 


f Rp + z s 
\ Rp + Zi 


2 


(63) 


or remain constant at the surface gravity, g s . The latter is 
computed by 


GM P 

(Rp + Zs) 2 


(64) 


where G is the gravitational constant, and z s is the altitude 
of the surface. 

We then compute the pressure as a function of altitude 
using 


( 65 ) 
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Table 2. Input and derived model atmosphere parameters for various planetary systems 



Earth 

temperate Jupiter 

HD189733b 

GJ1214b 

Parameters 

A 

B 

c 

D 

E 

F 

G 

H 

Input 

Composition a 

E 

E 

J 

J 

J 

J 

J 

J 

T (K) 

255 

255 

255 

255 

1090 

1090 

550 

550 

^Zatm (km) 

200 

200 

1000 

1000 

4500 

4500 

8800 

8800 

Az (km) 

0.01 

0.01 

0.01 

0.01 

0.05 

0.05 

0.05 

0.05 

Az samp (km) 

1 

1 

5 

5 

20 

20 

30 

30 

Az comp (km) 

0.01 

0.01 

0.05 

0.05 

0.05 

0.05 

0.1 

0.1 

Gravity ^ 

C 

V 

C 

V 

C 

V 

C 

V 

Derived 

H (km) 

7.4 

7.4 - 7.9 

35.7 

35.7 - 36.7 

174.5 

174.5 - 194.6 

239.7 

239.6 - 544.0 

Zs (km) 

-17.07 

-17.12 

-82.22 

-82.32 

-401.8 

-403.8 

-551.85 

-568.90 


a E: Earth-like; J: Jovian-like; See section 1.1.21 
^ V: varies with altitude; C: constant; See section 13.31 


starting with the surface at i = 0, and using the result 
of each iteration as input to the next, to the top of the 
atmosphere. The model atmosphere that is passed on to 
MAKEXOSHELL is sampled every Az samp , a much larger 
value than the thickness of a layer in MAKEXOSHELL’s 
computational grid, A z CO m P - Our choices of the vertical res¬ 
olution of the computational grid and of the model sampling 
grid are adapted to the scale height of the planetary atmo¬ 
sphere. Typically, on the order of 1000 computational layers 
fit inside one scale height, and the resulting model atmo¬ 
sphere is sampled roughly every 100 layers. 

We have chosen not to include the centrifugal force for 
several reasons. The most important one is that it requires 
knowledge of a planet’s rotation period, which for exoplanets 
is usually unknown. Coupled with the fact that the orienta¬ 
tion of the spin axis with respect to the observer is also un¬ 
known, one end up with several free parameters that can not 
be unambiguously determined. Since the effect of the cen¬ 
trifugal force is to reduce gravity preferentially at a planet’s 
equator, and hence increase an atmosphere’s scale height, 
that can not be told apart from a local increase in tem¬ 
perature or a decrease in the mean molecular mass of the 
atmosphere. Hence, considering the effects of the centrifugal 
force adds a level of complexity that is unwarranted given 
the aim of this paper. 

The isothermal temperature is chosen to be the mean 
planetary emission temperature, T e . This is given by 


Te 


1(1 — As) 



( 66 ) 


where the yet undefined parameters are the semi-major axis 
of the planetary orbit, a, the Bond albedo of the planet, As, 
and the effective temperature of the host star, T„. This emis¬ 
sion temperature is appropriate for a planet which thermally 
re-emits from its entire surface the absorbed stellar energy, 
and is often computed for various exoplanets over a breadth 
of possible Bond albedo, from 0 to 0.75 (see e.g. Charbon- 
neau et al. 2009 for GJ 1214b). In our simulations, we are 
choosing T e for a Bond albedo of 0.30, which is similar to 


the Bond albedo of Earth and the Jovian planets in our Solar 
system. 

Tabic [2] summarises the various input parameters of 
the model atmospheres considered in this paper. For each 
planet, we have one model where the gravity varies with 
altitude and one where it is constant, with otherwise iden¬ 
tical input parameters. All of the model atmospheres have 
a surface pressure of 10 atm (atmosphere), and a reference 
pressure of 1 atm. The atmospheric thicknesses have been 
chosen such that the top of the atmosphere where gravity 
varies with altitude have densities ranging between 10“ 12 
and 10“ 10 amagat. The table also lists the corresponding al¬ 
titude of the surface, and the range of density scale heights 
through the atmosphere. 


4 RESULTS AND DISCUSSIONS 

4.1 Validation of analytical expressions and ray 
tracing algorithm 

The column abundance from our analytical expressions, N a , 
is 

N a = N 0 + Vi + N 2 (67) 

expressed in equations l35l l37l and !38l and the analytical ray 
deflection, ui a , is 

UJa — CJO T (Ul -f- UJ 2 (68) 

expressed in equations 1391 1401 and 1411 The outputs from 
MAKEXOSHELL can include, or not, the contribution from 
the upper layers (See equations 1591 and 1601) . Comparison of 
one quantity, Q i, relative to another, Q 2 , is expressed by 
the percentage difference, 100(Qi — Q 2 )/Q 2 , either in terms 
of the absolute value or not. 

To verify the accuracy of our ray tracing algorithm, we 
compare the column abundance, N„, and ray deflection, w„, 
obtained with MAKEXOSHELL with those computed from 
our analytical expressions, for a planet with Earth’s radius 
and bulk mass (Rp = 6371.01 km; Mp = 5.9736 x 10 24 kg). 
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Figure 3. Percentage difference of the numerical, N n , relative to 
the analytical, N a , column abundance along a ray, as a function 
of the largest atmospheric density reached by that ray, for dif¬ 
ferent thicknesses (0.1, 0.01, and 0.002 km - see inset legend) of 
MAKEXOSHELL’s computational layers. That difference is also 
shown for a computational layer thickness of 0.01 km when the 
numerical contribution from the upper layers (see equation 1591 ) 
is not included (long-dashed line). The red curve shows our esti¬ 
mated contribution from the higher-order terms, N res t (see equa¬ 
tion 1451 . 


We use atmosphere model A, which obeys the same condi¬ 
tion as the assumptions used to derive our analytical expres¬ 
sions: isothermal temperature profile, constant composition 
with altitude, and constant gravity with altitude. The lower 
boundary (see section l2(2ll occurs at a density of 4.03 amagat 
(or pressure of 3.76 atm), at an altitude of —9.82 km. 

Figures [3] and [I] show the difference of the numerical out¬ 
put of MAKEXOSHELL relative to our analytical derivation 
for the ray column abundance as a function of the largest at¬ 
mospheric density reached by that ray. The nominal simula¬ 
tions use a thickness for the computational layers of 0.01 km 
and is shown by the solid line. For this vertical resolution, 
there is excellent agreement, better than 0.004 per cent, be¬ 
tween our numerical algorithm and our analytical expres¬ 
sions over a wide range of densities (from about 4 x 1CP 9 to 
6 x 10 -3 amagat). At these densities, the numerical column 
abundances are slightly lower than those from the analytical 
expressions, and the obtained accuracy seems to be limited 
by numerical instabilities. 

At lower densities, significant errors occur near the at¬ 
mospheric top boundary and the numerical algorithm finds 
larger column abundance than the analytical solution. This 
is due to our imperfect way of incorporating the contribu¬ 
tion of the upper layers. However, not including it causes 
much larger errors (see the long-dashed line). One solution 
is to adjust the top boundary to confine these errors to at¬ 
mospheric regions with densities low enough that they have 
negligible optical depths, and make no impact on the inter¬ 
pretation of remote sensing observations. 

At higher densities, the absolute difference in the col¬ 
umn abundance increases and seems to follow a density 
power law, until it reaches about 0.3 per cent around 1 am- 


Figure 4. Absolute percentage difference of the numerical, N n , 
relative to the analytical, N a , column abundance along a ray, as 
a function of the largest atmospheric density reached by that ray, 
for different thicknesses (0.1, 0.01, and 0.002 km - see inset legend) 
of MAKEXOSHELL’s computational layers. See the caption in 
Fig.0 



ABSOLUTE DEFLECTION DIFFERENCE (%) 


Figure 5. Absolute percentage difference of the numerical, LJ n , 
relative to the analytical, uj a , ray deflection, as a function of 
the largest atmospheric density reached by that ray, for differ¬ 
ent thicknesses (0.1, 0.01, and 0.002 km - see inset legend) of 
MAKEXOSHELL’s computational layers. That difference is also 
shown for a computational layer thickness of 0.01 km when the nu¬ 
merical contribution from the upper layers (see equation [60} is not 
included (long-dashed line). The red curve shows our estimated 
contribution from the higher-order terms, uj r est (see equation |46D . 


agat. In these density regions, the numerical results are still 
lower than the analytical expressions. We have investigated 
whether this could be due to an inadequate computational 
grid resolution, and have done calculations for vertical com¬ 
putational layer thicknesses, A z cornp , of 0.1 and 0.002 km 
(see the triple-dot dashed and the dot-dashed lines, respec- 
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tively, in Fig. [3] through [5]). The computational grid reso¬ 
lution does indeed play a role, as the differences between 
numerical and analytical results decrease markedly when 
changing the layer thickness from 0.1 to 0.002 km. Our 
choice of 0.01 km is a compromise between accuracy and 
computational memory requirements, with 740 to 790 com¬ 
putational layers per atmospheric scale height. 

For densities higher than about 1 amagat, the numer¬ 
ical results are higher than the analytical expressions be¬ 
cause the analytical expressions are missing the higher-order 
terms, and the difference increases dramatically with density 
as the lower boundary is approached. Although generally 
lower, it follows the same trend as our crude order of mag¬ 
nitude estimate of the contribution from these lrigher-order 
terms (see equation 1 45 II . shown by the red curve. 

Differences between the numerical results and the ana¬ 
lytical expressions of the ray deflection (see Fig. E]) are gener¬ 
ally higher than for column abundances, but follow roughly 
the same behaviour. The higher-order terms becomes non- 
negligible at densities starting around 0.15 amagat and the 
difference in ray deflection is about 1.6 per cent at 1 amagat. 

The excellent agreement between MAKEXOSHELL 
and our analytical expressions, in density regions where the 
higher-order terms are negligible, validates both methods si¬ 
multaneously. However, the analytical expressions can only 
be used to compute the column abundance along the ray, 
and its deflection with good accuracy up to densities of 
about 1 and 0.15 amagat, respectively. On the other hand, 
MAKEXOSHELL can be used in denser regions where it 
captures qualitatively the expected location of a singular¬ 
ity near the lower boundary, as well as under more general 
conditions such as when gravity changes with altitude and 
the temperature profile is not isothermal. Indeed, MAKEX- 
OSHELL’s only assumption is that the atmospheric density 
scale height is constant within each layer of the input model 
atmosphere, which is approximately true if the vertical sam¬ 
pling of the model atmosphere is high enough. 

Our analytical expressions are nevertheless very use¬ 
ful. They first provide a benchmark against which numeri¬ 
cal models can be compared. Indeed, we converged on the 
current numerical algorithm and vertical computational res¬ 
olution precisely because we could compare to our improved 
analytical expressions. They also allow derivation of analyt¬ 
ical expressions for other quantities which depend on the 
column abundance or the ray deflection, as well as allow 
propagation of input parameter uncertainties, both of which 
are crucial for atmospheric retrieval algorithms. 


4.2 Relevance to temperate exoplanets 

4-2.1 Column abundance and ray deflection 

With the validation of our analytical expression and our nu¬ 
merical model, we can ask to what degree does our analyti¬ 
cal expressions and numerical model differ from the simplest 
analytical expressions that are used in some exoplanet at¬ 
mosphere retrieval algorithm? The simplest analytical ex¬ 
pressions are noV^nHro for the column abundance, and 
vo \J2nro/H for the ray deflection. These are the leading 
term in our analytical expressions (see equations [43] and [44]). 
We consider this question both for the Earth-like planet that 


we have already described, as well as for a temperate Jovian- 
like planet. 

Our temperate Jovian-like planet combines atmo¬ 
spheric model C with Jupiter’s radius and bulk mass 
(Rp = 69911 km; Mp = 1.8986 x 10“' kg), and orbits 
a Sun-like star at Earth’s orbital distance. Since its density 
scale height is about five times that of our Earth-like planet, 
we can use a computational layer thickness of 0.05 km for 
MAKEXOSHELL and get similar agreements between the 
numerical results and our analytical expressions as for the 
Earth-like planet, because we have a comparable number of 
computational layers per scale height. 

Figures [6] and [7] show the percentage difference in the 
column abundance and the ray deflection, respectively, com¬ 
puted with our models, relative to those obtained with the 
simplest analytical expressions. This is shown both for our 
analytical expressions and our numerical algorithm, and for 
both type of planets considered. We are not showing the 
difference for densities lower than 10 -5 amagat because it 
is constant, except close to the top boundary where these 
quantities are overestimated by MAKEXOSHELL (such as 
in figures [3] through [5|. 

The curves for both types of planets are remarkably 
similar except at densities lower than 10~ 2 amagat. The 
difference in those density regions comes from the Co and 
the Do coefficients in equations 1431 and 1441 which are power 
series of the ratio ( H/r ). The similarities of these curves for 
both planets at deeper densities is fortuitous as the factor of 
2 difference in ( r/H ) is almost compensated by differences 
in refractivities, so that the dimensionless parameter ( rv/H ) 
of these two planets is roughly the same at similar densities. 
Consequently, the lower boundary occurs on the temperate 
Jupiter at a density of 4.15 amagat (or pressure of 3.87 atm), 
only marginally denser than on the Earth-like planet, at an 
altitude of —52.32 km. 

For both planets, the simplest analytical expressions 
underestimate the column abundance and ray deflection by 
about 1 per cent around densities of 0.1 amagat («0.1 atm), 
10 per cent around 1 amagat (~1 atm), and this error in¬ 
creases dramatically as rays approach the lower boundary, 
where (rv/H) is equal to one. Thus, in the denser atmo¬ 
spheric regions, the simple expressions are quite inaccurate. 

4-2.2 Effective exoplanet radius 

Can these regions be observed in a transiting exoplanet, or 
are the optical depths there so high that these differences 
have no impact on an exoplanet’s effective radius? When 
stellar limb darkening is ignored, the effective radius, R e ff, 
of an exoplanet occulting the centre of its host star is given 
by 

Riff = b 2 max - 2^ e~ T bdb (69) 

where r = aN is the optical thickness traversed by a ray, 
a is an average atmospheric extinction cross section, and 
bmax is the impact parameter of rays grazing the top of the 
atmosphere. 6 m i n is the impact parameter of rays grazing ei¬ 
ther the critical altitude or the planetary surface, whichever 
is larger. For a cloud-free homogeneous isothermal atmo¬ 
sphere, when one ignores refraction, the effective radius of 
an exoplanet occurs where grazing rays have traversed an 
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Figure 6. Percentage difference of the numerical (solid line) and 
analytical (dashed line) ray column abundance relative to that 
from the analytical expression nQy/2nHro , as a function of the 
largest density reached by a ray, for Earth-like (black) and tem¬ 
perate Jovian-like (red, see text) planets. 

optical thickness of about 0.56 (Lecavelier des Etangs et al. 
2008). This result is sens itive roughly to four atmospheric 
scale height (|Gri ffithl20l4l . where the atmospheric transmis¬ 
sion e~ T changes from 0.05 to 0.95, and column abundances 
are 5.4 and 0.1 times larger than at the effective radius, 
respectively. 

As opacities change with wavelength, so does the ef¬ 
fective radius of the exoplanet. A clear Rayleigh scattering 
atmosphere free of other sources of opacity produces the 
smallest effective radius. The average Rayleigh cross sec¬ 
tion are about 4.0 x 10 -28 cm 2 and 7.5 x 10~ 29 cm 2 at 
1 pm, for an Earth-like and Jovian-like composition, re¬ 
spectively. The smallest effective radius at that wavelength 
then occurs where the column abundance along the ray are 
1.4x 10 2 ' cm -2 and 7.4x 10 27 cm -2 , respectively. The largest 
effective radius of a cloud-free Earth in the NIR occurs in 
the core of a CO 2 band at 4.3 pm at a column a bu ndanc e 
of about 3.5 x 10 23 cm' 2 dBetremieux fe Kalteneggerll2014l) . 
which we will use as a rough guideline of the largest value in 
the NIR of the effective radius of a planet for the discussion 
throughout this paper. If one is interested in far ultraviolet 
spectral features, higher altitude regions of the atmosphere 
are actually of importance. 

It is important to stress that the various column abun¬ 
dances quoted above are meant only as indicators to roughly 
translate the size of spectral features produced by Earth’s 
atmosphere into equivalent features in other types of atmo¬ 
sphere, as well as get a feel for the range of atmospheric den¬ 
sities that contribute significantly to the creation of spectral 
features. They are in no way a substitute for actual radia¬ 
tive transfer calculations. In figures ISlthrough l 151 the triple- 
dot-dashed lines bracket the grazing altitude corresponding 
to the smallest and largest effective radius using the sim¬ 
ple column abundance criterion described above, as well as 
the range of ray deflection in that region. The thick blue 
vertical lines bracket the range of column abundances to 
which these spectral features are sensitive, from about 10 


Figure 7. Absolute percentage difference of the numerical (solid 
line) and analytical (dashed line) ray deflection relative to that 
from the analytical expression , as a function of the 

largest density reached by a ray, for Earth-like (black) and tem¬ 
perate Jovian-like (red, see text) planets. Note that the two min¬ 
ima around 10 —3 amagat indicate a transition to negative differ¬ 
ences which occur at lower densities. 

times lower than at the largest effective radius to about 5.4 
times larger than at the smallest effective radius. If our re¬ 
fractive model modifies substantially the column abundance 
within this region, then the cond itions for the validity of 
iLecavelier des Etangs et al.l (120081 ) ’s result breaks down and 
the effective radius of the exoplanet will hence be modified. 

Figures [8] through |TT] show the ray column abun¬ 
dance and deflection computed with MAKEXOSHELL for 
our isothermal Earth-like and temperate Jovian-like atmo¬ 
sphere, as a function of a ray’s grazing altitude. These are 
shown both for atmospheric models B and D where gravity 
changes with altitude (solid line) and for models A and C 
where it is kept constant at its surface value (dashed line). 
Note that in all these simulations, the contribution from the 
upper layers are included. 

For both of these temperate planets, the density scale 
height is much smaller than their radius, so that the varia¬ 
tion of gravity with altitude does not make much of a change 
to the column abundances or ray deflection in the atmo¬ 
spheric regions that create spectral features. However, these 
regions include high density regions where our treatment of 
refraction computes column abundances which are at least 
4 times larger than thos e co mp uted wit h the simple analyt¬ 
ical expressions used by ILecavelier des Etangs et al] (2003) • 
Hence, our treatment of refraction will modify the value of 
a transiting temperate exoplanet’s effective radius at wave¬ 
lengths where opacities are low. 

4-2.3 Effective versus actual scale height 

Furthermore, the effective radius depends on the variation 
of the column abundance along the ray with respect to its 
impact parameter, and not to the grazing altitude of a ray. 
When plotted versus the impact parameter (or in this case 
versus b — Rp, which is equivalent), both the column abun- 
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Figure 8. Ray column abundance as a function of the grazing al¬ 
titude of a ray, for the Earth-like isothermal atmosphere described 
in sections 14.11 and 14.21 when gravity varies with altitude (solid 
line), and when gravity is constant (dashed line). The red line is 
identical to the solid line except that it is plotted as a function 
of the projected distance to the reference planetary radius with 
respect to the observer. The triple-dot-dashed lines bracket the 
likely effective radii of spectral features, and the vertical blue line 
bracket the range of column abundances to which these spectral 
features are sensitive (see section f4.2.211 . 



Figure 9. Ray deflection as a function of the grazing altitude 
of a ray, for the Earth-like isothermal atmosphere described in 
sections 1 1.1 1 and 14.21 See the caption in Fig. [8] 

dance and the ray deflection display a behaviour in the 
denser regions which is characteristic of a much smaller scale 
height (as shown by the solid red curve). The effective ra¬ 
dius of the exoplanet reaches asymptotically a value, and 
thus effectively mimics a surface, which is characteristic of 
a much smaller density than the highest density that can 
be probed. For our isothermal Earth atmosphere, the lower 
boundary, located at a density of 4.03 amagat, maps to an 
impact parameter 2.4 km below our 1 atm reference. If one 



Figure 10. Ray column abundance as a function of the grazing 
altitude of a ray, for the temperate Jovian-like isothermal atmo¬ 
sphere described in section 11.21 See the caption in Fig. [8] 



Figure 11. Ray deflection as a function of the grazing altitude 
of a ray, for the temperate Jovian-like isothermal atmosphere de¬ 
scribed in section T4.21 See the caption in Fig. [8] 

ignores this effect, one concludes that this impact parame¬ 
ter correspond to a density of only about 1.38 amagat. In 
the temperate Jovian-like atmosphere, the lower boundary, 
located at a density of 4.15 amagat maps to an impact pa¬ 
rameter 12.7 km below our 1 atm reference, or an apparent 
density of 1.43 amagat. 

These effects are due to the refractive deflection of the 
ray, and an altitude interval maps to increasingly smaller 
impact parameter interval for larger densities. Indeed this 
can be seen mathematically by differentiating equation [T| 
with respect to ro, and then using equations 0 evaluated at 
ro to obtain, 

k- ,+ "(*-3)- (™) 

A ray that traverses the atmosphere must satisfy equation [9] 
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Table 3. Given and derived parameters for HD 189733b 




References 

a 

Parameters 

B2005 

T2008 

S2010 

Given 

R* (R©) 

0.760 

0.756 

0.752 

M, (M©) 

0.820 

0.806 

0.840 

T, (K) 

5050 

5040 

5050 

R p (R j) 

1.26 

1.138 

1.151 

M P (Mj) 

1.15 

1.144 

1.150 

a (AU) 

0.0313 

0.0310 

0.03142 

Derived 

p (x 10 3 kg m 3 ) 

0.763 

1.030 

1.000 

g s (m s“ 2 ) 

18.76 

22.90 

22.50 

T e (K), A s = 0 

1200 

1200 

1191 

T e (K), As = 0.3 

1098 

1098 

1090 

T e (K), As = 0.5 

1009 

1009 

1002 

T e (K), As = 0.75 

848.6 

848.8 

842.5 

UJ C (°) 

7.572 

7.504 

7.381 


a B2005:lBouchv et all J2005h : T2008: lTorres et al.l [|2008h : S2010: 
ISouthworthl d201 (1 ) 

which results in ( db/dro) > 0. In the upper regions of the at¬ 
mosphere where refractivities tend to 0, this derivative tends 
to 1, and altitude intervals are almost identical to impact pa¬ 
rameter intervals. Near the lower boundary, this ratio tends 
to 0, so that an altitude interval maps to a vanishingly small 
impact parameter interval, thus forming what amounts to a 
refractive boundary layer. 

For our isothermal Earth atmosphere, the effective scale 
height of the observed column abundance has an average 
value of 4.3 km between 1.1 and 1.5 amagat, 1.8 km be¬ 
tween 2.1 and 2.9 amagat, and only 0.32 km between 2.9 
and 4.0 amagat, in contrast to the actual scale height of 
the atmosphere of about 7.4 km. A similar rapid decrease in 
the effective scale height occurs for the temperate Jovian at¬ 
mosphere. When expressed in term of the planet’s effective 
radius, the size of spectral features in transmission spec¬ 
troscopy scale to first order with the effective scale height. 
This implies that features that originate in those dense re¬ 
gion, such as those from dimers or collision induced absorp¬ 
tion, will appear that much smaller and will be harder to 
detect. 

How much of an impact this effect has on the transit 
lightcurve depends on the part of the transit. When the 
planet occults the centre of its star, one can not observe 
atmosph eric regions d enser than t he cri tical density of the 
planet dBetremieux fe Kalteneggerll2014l) . During the course 
of the transit, the symmetry is broken, and one part of the 
planet’s limb can be probed to higher densities while the 
opposite limb can only be probed to lower densities. If the 
critical density of the planet occurs at a density where the 
effective scale height is not significantly different from the 
actual scale height, then this apparent decrease of the scale 
height of the denser regions may not be observed, except 
possibly at the edges of the transit. 

For a given column abundance, the Earth-like planet 
bends radiation more than the temperate Jovian-like planet. 


For instance, at a column abundance of 3.5 x 10 23 cm~ 2 , 
the ray deflection are about 2.8 x 10 -4 and 2.5 x 10 -5 de¬ 
gree, respectively, roughly a factor of 11 difference. The crit¬ 
ical density to which one can probe the atmosphere during 
an exoplanet transit is higher for the temperate Jovian-like 
planet. Around a G2 star, the critical densities are 0.22 and 
0.36 amagat for the Earth-like and the temperate Jovian-like 
planets, respectively. Around an M9 star (R* = 0.08 R©; 
T* = 2300 K), the planet must have a semi-major axis of 
0.0127 AU, in order to receive the same flux as a planet 
1 AU from a G2 star. Since the planetary radii are compa¬ 
rable to the stellar radius, the critical deflections of the two 
planets are markedly different (1.87°and 3.79°, respectively) 
and so are the critical densities (1.32 and 3.00 amagat, re¬ 
spectively), which are high enough that the effective scale 
height of the atmosphere is significantly different from the 
actual scale height. 

As our treatment of refraction makes a difference for 
temperate planets, we will now investigate whether this is 
also the case for hotter exoplanets which are more eas¬ 
ily observed. As an example, we will consider two planets 
that have already been repeatedly observed: the hot Jupiter 
HD189733b and the super-Earth GJ1214b. 

4.3 HD189733b 

The NASA Exoplanet Archive dAkeson et alj l2013li lists 
three references for the system parameters o f HD18 9733b 
and its host stars: its discovery bv lBouch v et al. l (120051). an d 
furt her observations an d reanalvsis bv lTorres et al.l (12008 J. 
and ISouthworthl (|2010| J. Table [3] lists some of those param¬ 
eters, as well as the mean mass density, p, surface gravity, 
g s , the mean emission temperature for various Bond albedo, 
T e (As), of the planet, and the critical deflection, w c , of the 
planet-star system (see Betremieux & Kaltenegger 2014 for 
its definition), that we derive from these parameters. Note 
that we list the temperature for albedos between 0 and 0.75, 
the same range conside r ed bylCharbonneau et al.l d2009lj and 
lAnglada-Escude et al.l (120131 ) for GJ1214b. Since its discov¬ 
ery, the planetary and stellar parameters of the HD 189733 
system have evolved quite substantially, but nevertheless 
the mean mass densit y and surface grav ity o f the planet 
from t he parameters of lTorres et al.l (I2008I J and ISouthworthl 
d20ld) are very similar. What is very uncertain remains the 
temperature of the atmosphere, as this depends on the com¬ 
position of minor species and whether clouds and aerosols 
are present. Hence, it is not worth to get into a discussion 
about whose planetary and stellar parameters are the best, 
as they do not dominate the uncertainties of the atmospheric 
properties, and we will merely use the revised parameters of 
ISouthworthl (l2010l ) to get a sense of what effects are of im¬ 
portance. 

We combine the planetary parameters (radius and 
mass) of HD189733b listed under S2010 in Tablc[3] with at¬ 
mosphere models E and F. Figures fl2l and fl3l show the ray 
column abundance and deflection as a function of the graz¬ 
ing altitude of a ray. Altitude regions where the effective 
scale height decreases dramatically is below the modelled 
regions. Just above our 10 atm “surface”, the average effec¬ 
tive scale height between 1.82 and 2.50 amagat is 146.1 km 
while the actual density scale height is 174.5 km. The den¬ 
sities at a given pressure are about a factor of 4.3 lower 
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Figure 12. Ray column abundance as a function of the graz¬ 
ing altitude of a ray for the Jovian-like isothermal atmosphere of 
HD189733b described in section f4. 3 1 using the S2010 parameters 
of Table [3] See the caption in Fig. [8] 


Figure 13. Ray deflection as a function of the grazing altitude 
of a ray for the Jovian-like isothermal atmosphere of HD189733b 
described in section T4.31 using the S2010 parameters of Table [3] 
See the caption in Fig. \8\ 


than for our temperate planets because the temperature is 
so much higher. Consequently, (rv/H) never reaches unity, 
even at the “surface” where it is only about 0.14, and col¬ 
umn abundance and ray deflection are increased by about 
6 per cent by our treatment of refraction compared to the 
simplest analytical expressions. 

The critical deflection of the planet-star system is about 
7.4° (see bottom of Table 0 , larger than the surface deflec¬ 
tion of about 1.0°. Consequently, the critical altitude will 
occur below our “surface”, and thus below the regions to 
which spectral features are sensitive (see blue vertical lines 
in Fig. ED. Hence, refraction does not prevent any part of 
the atmosphere from being observed at wavelengths below 
1 pm. Although Rayleigh scattering cross sections decrease 
with the fourth power of the wavelength, many molecular 
absorption lines exist in the infrared, and are usually the 
dominant source of opacity. 

Given our simple column abundance criterion (see sec¬ 
tion [472T2J , the larger effective radius of the non-constant 
gravity model is 60 km larger than the constant gravity 
model, or about a third of a scale height, which translates 
into an extra stellar flux drop of 3.5 x 10 -3 per cent during 
transit. Since this is smaller than uncertainties in observa¬ 
tions of HD189733b to date (see e. g. Fig. 4 in McCullough et 
al. 2014), this effect is not important for current exoplanet 
observations, but might be for future James Webb Space 
Telescope observations. 

4.4 GJ1214b 

The NASA exoplanet archives list a few references for the 
system pa rameters of GJ1214b and its host star: its dis¬ 
covery by Charbonncau et_alj ( 200911, and fu r ther obser¬ 
vations by Carter et alJ (1201 if) . HarpsOe et al.l (120131 ). and 
lAnglada-Escude et al.l (120131 ). Table [4] lists the given pa¬ 
rameters from these different sources, as well as the mean 
mass density, surface gravity, average emission tempera¬ 
ture, and critical deflection, that we derive from them (see 


section m for sy mbol defin itions). Of these observations, 
lAnglada-Escude et al.l (120131 ) uses recent parallax observa¬ 
tions to re-evaluate the distance of the system, and com¬ 
bines it with new radial velocity measurements and transit 
observations to refine the system parameters. The derived 
planetary parameters are very s imilar to those previously ob¬ 
tained bv iHarpsOe et all (120131 ). namely a low mean density 
and small surface gravity, but with higher temperatures than 
previously found. In this section, we will use their “maxi¬ 
mum probability” planetary parameters (radius and mass) 
listed under AE2013 in Table [4] with atmospheric models G 
and H. 

Figures [14] and [15] show the ray column abundance and 
deflection as a function of the grazing altitude of a ray. At 
the “surface”, ( rv/H ) is 0.044, smaller than the value for 
HD189733b, and our treatment of refraction makes even less 
of a difference than for HD189733b. Indeed, column abun¬ 
dances at the “surface” are increased by about 3 per cent, 
and ray deflections by 1 per cent, by our treatment of re¬ 
fraction. The ray deflection at the surface is about 0.76°, 
far less than the critical deflection of 4.4°, and the critical 
altitude is below our “surface”. Hence, for a Jovian composi¬ 
tion, refraction does not prevent any part of the atmosphere 
of GJ1214b from being observed. 

However, the combination of a high temperature with a 
low gravity and a low mean molecular mass creates a scale 
height which is a non-negligible fraction of the planetary 
radius. As a result, the variation of gravity with altitude 
makes quite a difference to the density altitude profile and 
to the resulting column abundances. Indeed, in order for 
atmospheric model H to have comparable densities at the 
upper boundary as models B, D, and F, its atmospheric 
thickness must be significantly larger (see Tabled - Ignoring 
the variation of gravity with altitude causes a drop in the 
column abundance of 5 orders of magnitude at that upper 
boundary, whereas it is less than 1 order of magnitude for 
the other models. 

The size of spectral features increases by about 655 km 
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Table 4. Given and derived parameters for GJ1214b 


Parameters 

Ch2009 

Ca2011 
Method A 

References a 
Ca2011 H2013 

Method B 

AE2013 

Maximum 

probability 

AE2013 

Expected 

value 

Given 

R* (R©) 

0.211 

0.210 

0.179 

0.216 

0.213 

0.211 

M* (Mg) 

0.157 

0.157 

0.156 

0.150 

0.176 

0.176 

T* (K) 

3026 

3170 

3170 

3026 

3252 

3252 

Rp (Re) 

2.678 

2.65 

2.27 

2.85 

2.80 

2.72 

M P (M s ) 

6.55 

6.45 

6.43 

6.26 

6.26 

6.19 

a (AU) 

0.01439 

0.01437 

0.01225 

0.01411 

0.01449 

0.01435 

Derived 

p (x 10 3 kg m 3 ) 

1.88 

1.91 

3.03 

1.49 

1.57 

1.70 

g s (m s“ 2 ) 

8.97 

9.02 

12.25 

7.57 

7.84 

8.22 

T e (K), A s =0 

559 

584 

584 

571 

601 

601 

T e (K), A s =0.3 

511 

535 

535 

522 

550 

550 

T e (K), As = 0.5 

470 

491 

491 

480 

506 

506 

Te (K), As =0.75 

395 

413 

413 

404 

425 

425 

Ulc (°) 

4.367 

4.349 

4.351 

4.578 

4.394 

4.386 


a Ch2009: ICharbonneau et al.1 120091); Ca2011: I Carter et al.l d201lh : H2013: lHarpsge et al.l 
d2013h ; AE2013: lAnglad^Escude et all ll2013h 




Figure 14. Ray column abundance as a function of the graz¬ 
ing altitude of a ray for the Jovian-like isothermal atmosphere of 
GJ1214b described in section 14.41 using the AE2013 maximum 
probability parameters of Table [4] See the caption in Fig. [8] 


Figure 15. Ray deflection as a function of the grazing altitude 
of a ray for the Jovian-like isothermal atmosphere of GJ1214b 
described in section [A4] using the AE2013 maximum probability 
parameters of Table [T] See the caption in Fig. [8] 


when considering the effect of a non-constant gravity with 
altitude, or 0.44 per cent of the stellar radius, which is com¬ 
parable to uncertainties in recent observations (see e.g. Fig. 4 
in de Mooij et al. 2013, and references therein). However, one 
should then also expect spectral features to change the effec¬ 
tive radius of the planet by about 3050 km, or about 2 per 
cent of the stellar radius, rather than the 0.6 per cent that 
has been observed to date. 

We will not try here with our illustrative model to solve 
the mystery of GJ1214b’s atmospheric composition but we 


will merely point out, as done by many previous analyses 
(see e.g. Bean et al. 2010; Croll et al. 2011; Crossfield, Bar¬ 
man & Hansen (2011); Frame et al. 2013; de Mooij et al. 
2013; Kreidberg et al. 2014), that a haze or cloud layer 
is a possible solution. Indeed, the top of that layer would 
probably determine the lowest altitude that could be ob¬ 
served, rather than Rayleigh scattering, decreasing the alti¬ 
tude range available for formation of spectral features. In- 
cidently, our zero altitude reference pressure would also be 
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correspondingly lower, shifting the densities in our atmo¬ 
spheric model to smaller heights. 


5 SUMMARY AND CONCLUSIONS 

We have derived new analytical expressions for the column 
abundance along a ray refracted by a planetary atmosphere 
with a constant density scale height, as well as the deflection 
experienced by that ray. These expressions were obtained by 
doing a Taylor power-series expansion of the exact integral 
formulation to the second order of the highest refractivity 
reached by the refracted ray. Surprisingly, the resulting an¬ 
alytical expression is not dependent on powers of the refrac¬ 
tivity, but rather on powers of the dimensionless parame¬ 
ter ( rv/H ). Although the refractivity is typically small, this 
dimensionless parameter is not necessarily so, and higher- 
order terms are dominant in a temperate atmosphere when 
rays reaches densities larger than about 1 and 0.15 amagat 
for the column abundance and the deflection, respectively. 
We have thus also derived an analytical expression that gives 
a rough estimate of the contribution from these higher-order 
terms. 

Since the scale height of most atmospheres are usually 
not constant, we have also developed a new numerical algo¬ 
rithm, MAKEXOSHELL, to trace rays through a spherically 
symmetric atmosphere given the STP refractivity of the at¬ 
mosphere and an arbitrary one-dimensional temperature- 
pressure profile. We have compared the results of numeri¬ 
cal simulations to values obtained with our analytical ex¬ 
pressions for a 10 atm 255 K isothermal Earth-like atmo¬ 
sphere on a planet with a size and mass identical to Earth’s. 
We have found that there is a better than 0.004 per cent 
agreement for densities up to 6 x 10~ 3 amagat, both for the 
computed column abundance and the ray deflection, except 
near our upper atmospheric boundary where densities are 
extremely low, and in the denser regions where the higher- 
order terms dominate. The excellent agreement, between our 
analytical expressions and numerical results, validates both 
methods simultaneously. 

We have built a few simple isothermal atmospheric 
models (temperate Earth-like planet, temperate Jovian-like 
planet, GJ1214b, and HD189733b) to determine the type of 
transiting exoplanets where our treatment of refraction, as 
well as the variation of gravity with altitude, makes an im¬ 
pact on their effective radius. The effects of a non-constant 
gravity with altitude are more pronounced for atmospheres 
where the atmospheric scale height is not much smaller than 
the planetary radius. The combination of a hot tempera¬ 
ture with a low gravity and mean molecular mass, creates 
the proper condition in a Jovian-composition GJ1214b. If 
GJ1214b’s atmosphere had a cloud or haze-free Jovian com¬ 
position, the difference in density profile due to the varia¬ 
tion in gravity is sufficiently large that the resulting change 
in effective radius at high opacities is comparable to cur¬ 
rent observational uncertainties. It is far less important for 
HD 189733b, and negligible for both temperate planets. 

The difference in the computed column abundance and 
the ray deflection, between our numerical treatment of re¬ 
fraction and simple analytical expressions used in current 
exoplanet atmosphere retrieval algorithms, increases with 
the dimensionless parameter (rv/H), and thus with the 


largest atmospheric density reached by a ray. We find that 
for hot Jupiters, such as HD189733b, as well as a Jovian- 
composition GJ1214b, atmospheric regions to which spec¬ 
tral features are sensitive are not dense enough for our im¬ 
proved treatment of refraction to make much of an impact. 
However, for temperate isothermal atmospheres of Earth¬ 
like and Jovian-like planets, the difference is about 1 per 
cent at 0.1 amagat, 10 per cent at 1 amagat, and increase 
dramatically at larger densities, which can impact spectral 
features. 

The column abundance and ray deflection tend toward 
infinite values, thus forming a refractive boundary layer, 
as the deepest region that can be probed by transmission 
spectroscopy and stellar occultation is approached. This 
lower boundary occurs where the radius of curvature of a 
grazing ray exactly matches the radial position of the ray 
with respect to the centre of the planet, and is only depen¬ 
dent on the density profile of the planet’s atmosphere. Rays 
that reach this boundary will spiral deeper into the atmo- 
sp here and b e absorbe d. The cri tical altitude discussed in 
iBetremieux fe Kalteneggeil d2014ll . can never be lower than 
this lower boundary no matter how close the exoplanet is to 
its star. 

Combined with the bending property of the atmo¬ 
sphere, temperate atmospheric layers with densities larger 
than 1 amagat have an apparent scale height significantly 
smaller than the actual density scale height of the atmo¬ 
sphere. This effective scale height tends to zero as rays ap¬ 
proach the lower boundary. In transmission spectroscopy, 
spectral features scale to first order with scale height when 
expressed in term of the effective planetary radius. Hence, 
spectral features that originate in these dense regions, such 
as those from dimers and collision-induced absorption, will 
be much harder to detect. Furthermore, this interesting ef¬ 
fect effectively mimics a surface so that even temperate gi¬ 
ant planets seem to have a surface in spectral regions where 
opacities are low. This effect can only be seen for exoplanets 
that are close enough to their star that the critical altitude 
of the planet-star system approaches this lower boundary. 
Given that it is easier to detect a planet by transit spec¬ 
troscopy around cooler stars, and that as a result many ex¬ 
oplanet searches now specifically target M-dwarf stars, this 
is more likely to be indeed the case. 
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